%%This file aims to create IVs---averaging market share and market dominance
%indicators across all outside associated markets for the longer panel


clear
tic

textFileName1=sprintf('sysyr.raw'); %the input file: the [hsa year sysid] for involved multi-region system
textFileName2=sprintf('input4createIV_oldyrs.raw'); %the input file of basic info of all hospitals
sysyr=dlmread(textFileName1);  %matrix: [hsa year sysid]
H=dlmread(textFileName2);  %matrix: [newid year prechs adj_vnd hsa sysid bdtot]
choice=1:13;
ns=numel(choice);


%%%%%%%%%%%%%%%%%%%%
%The first way is to find out the largest HSA according to how many
%member hospitals locate in.
exo_mktdom=zeros(size(sysyr,1),ns);
exo_mktsh=zeros(size(sysyr,1),ns);
for i=1:size(sysyr,1)
    tempsys=H(H(:,6)==sysyr(i,3)&H(:,2)==sysyr(i,2),:);
    %%%%%%%%%%
    %The following is to compute the exogenous big hsa for a hospital
    %system
    syshsa=sort(tempsys(:,5));
    unique_syshsa=unique(syshsa);
    occur=histc(syshsa,unique_syshsa);
    bighsa=unique_syshsa(occur==max(occur));
    exo_bighsa=bighsa;
    exo_bighsa(exo_bighsa(:)==sysyr(i,1))=[];
    %if min(abs(exo_bighsa-sysyr(i,1)))==0
     %   exo_bighsa=[];
    %end
  if isempty(exo_bighsa)==0 
      part_mktdom=zeros(numel(exo_bighsa),ns);
      part_mktsh=zeros(numel(exo_bighsa),ns);
      for m=1:numel(exo_bighsa)
          prechs=H(H(:,5)==exo_bighsa(m)&H(:,2)==sysyr(i,2),3);
          prechs(prechs==1)=[];
          if isempty(prechs)==0
          occur=histc(prechs,choice);
          occur=reshape(occur,1,[]);           
          part_mktsh(m,:)=occur/sum(occur); %a vector 1 by 13
          temp_prechs=prechs;
          temp_prechs(prechs==2|prechs==13)=[];
          if isempty(temp_prechs)==0
             occur=histc(temp_prechs,choice)';
             mktdom=choice(occur==max(occur));
             Tchoice=choice';                     
             vmktff=max((~(Tchoice(:,ones(size(mktdom,2),1))-mktdom(ones(size(choice,2),1),:))),[],2);
             part_mktdom(m,:)=vmktff';            
          end 
          end
      end
      exo_mktdom(i,:)=mean(part_mktdom,1);
      exo_mktsh(i,:)=mean(part_mktsh,1);
  end 
end


iv_mktdom=[sysyr(:,1:3),exo_mktdom];
iv_mktdom(~any(iv_mktdom(:,4:end),2),:)=[];  %### UPDATED on 07/25/2019
iv_mktsh=[sysyr(:,1:3),exo_mktsh];
iv_mktsh(~any(iv_mktsh(:,4:end),2),:)=[];  %### UPDATED on 07/25/2019

iv_mktdom=sortrows(iv_mktdom,[1,2,3]);
iv_mktsh=sortrows(iv_mktsh,[1,2,3]);

hsayr=unique(sortrows(iv_mktdom(:,1:2),[1,2]),'rows');

output_ivdom=zeros(size(hsayr,1),ns);
output_ivmsh=zeros(size(hsayr,1),ns);
for i=1:size(hsayr,1)
    output_ivdom(i,:)=mean(iv_mktdom(iv_mktdom(:,1)==hsayr(i,1)&iv_mktdom(:,2)==hsayr(i,2),4:end),1);
    output_ivmsh(i,:)=mean(iv_mktsh(iv_mktsh(:,1)==hsayr(i,1)&iv_mktsh(:,2)==hsayr(i,2),4:end),1);   
end


textFileName11=sprintf('iv_dom.csv');
dlmwrite(textFileName11,[output_ivdom,hsayr],'delimiter', ',', 'precision', 8)
textFileName12=sprintf('iv_msh.csv');
dlmwrite(textFileName12,[output_ivmsh,hsayr],'delimiter', ',', 'precision', 8)







